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INTRODUCTION 


The  flow  around  airfoils,  struts  and  obstacles  placed  in  an  approaching  stream 
with  a  non-uniform  velocity  is  a  fundamental  three-dimensional  viscous  I  low  of 
considerable  importance.  Examples  of  this  type  of  flow  include  the  flow;;  near  an 
aircraft  wing/fuselage  junction  and  near  a  submarine  hull /sail  junction.  Another 
example  occurs  in  axial  compressors  and  turbines,  where  the  houndarv  lovers  which 
develop  on  the  annular  surfaces  of  the  axial  flow  passage  encounter  rows  of  stat ionarv 
and  rotating  blades.  Other  examples  include  the  flow  of  a  river  past  a  bridge  pier 
or  other  underwater  structure  and  other  flows  past  surfaces  where  a  protuberance  is  of 
necessity  present.  Examples  can  even  be  found  in  such  exotic  devices  as  fluidic 
control  systems. 

The  feature  common  to  all  such  flows  is  that  a  non-uniform  velocity  in  an  approach¬ 
ing  boundary  layer  meets  a  local  region  of  adverse  pressure  gradient  due  to  the  blockage 
effect  of  the  obstruction.  This  produces  a  complex  flow  field  consisting  of  houndarv 
layer  separation  and  the  formation  of  one  or  more  horseshoe  vort ices  around  the 
obstruction.  Regions  of  interest  include  both  the  horseshoe  vortex  region  at  the 
lending  edge  and  the  corner  flow  region  downstream  of  the  leading  edge.  In  the  case  of 
an  airfoil  or  strut,  this  corner  flow  region  downstream  of  the  leading  edge  contains 
streamwise  vortices  whicli  affect  both  the  performance  of  the  airfoil  or  strut  them¬ 
selves  and  also  the  performance  of  other  flow  devices  located  downstream.  Tn  flows 
where  heat  transfer  is  critical,  the  presence  of  the  streamwise  vortices  causes  a 
local  thinning  of  boundary  layers  and  transport  of  f roe  stream  fluid  toward  the  wall 
in  portions  of  the  flow  region.  This  in  turn  serves  to  increase  local  heat  transfer 
rates.  Order-of-magnitude  increases  and  local  heat  transfer  rates  have  been  observed 
experimentally. 

Although  the  horseshoe  vortex  flow  in  the  region  of  the  leading  edge  is  of 
interest  in  its  own  right,  it  is  also  of  considerable  importance  in  connection  with 
the  corner  flows  which  occur  downstream.  Corner  flows  are  often  viewed  as  an  isolated 
problem  without  giving  much  attention  to  the  leading  edge  region  where  the  corner 
flows  originate.  In  cases  where  leading  edge  effects  are  minimal  or  localized 
(e.g.  sharp  or  cusped  leading  edges  placed  at  zero  effective  incidence  to  the  approach 
flow),  this  approach  may  be  sufficient.  It  seems  doubtful,  however,  whether  corner 
flows  in  general  can  be  divorced  from  the  flow  in  the  leading  edge  region  where  the 
corner  tlows  are  formed,  if  for  no  other  reason  than  to  determine  appropriate  initial 
conditions  for  a  corner  flow  analysis. 


Several  experimental  flow  visualization  studies  of  the  horseshoe  vortex  flow 
have  been  performed,  and  these  have  established  that  the  flow  consists  of  a  three- 
dimensional  boundary  layer  separation  in  front'  of  the  obstruction  followed  by  a  vortex 
flow  which  wraps  around  the  obstruction.  Little  is  available  in  the  way  of  detailed 
flow  measurements  (including  all  three  velocity  components)  for  the  horseshoe  vortex 
problem,  particularly  downstream  of  separation.  The  only  available  measurements  which 
include  the  secondary  velocity  components  appear  to  be  the  recent  measurements  of  Mehta, 
Shabaka  and  Bradshaw  [1]  for  an  unswept  leading  edge  at  zero  incidence,  taken  in  the 
corner  flow  region  well  downstream  of  the  leading  edge,  where  the  secondary  flows  are 
relatively  weak.  Although  these  measurements  do  not  cover  the  leading  edge  region  where 
the  horseshoe  vortex  is  formed,  the  measurements  are  nonetheless  helpful  in  assessing 
computational  flow  predictions. 

Previous  analytical  studies  have  considered  the  three-dimensional  boundary  layer 
flow  upstream  of  separation  (e.g.,  Dwyer  [2])  and  have  used  rotational  inviscid  flow 
theory  to  estimate  secondary  flows  (Hawthorne  [3]).  The  authors  [4]  have  computed 
solutions  of  the  compressible  Navier-Stokes  equations  for  laminar  flow  near  an  elliptical 
leading  edge,  including  the  leading  edge  flow  separation,  horseshoe  vortex  region,  and 
the  formation  of  corner  flows.  Solutions  are  given  in  [4]  for  both  zero  and  five 
degree  angles  of  incidence.  In  the  present  study,  this  approach  is  extended  to 
encompass  turbulent  flow  and  more  general  geometries.  Solutions  are  given  for  tur¬ 
bulent  horseshoe  vortex  flow  past  both  unswept  and  45  degree  swept  elliptical  leading 
edges  affixed  to  a  fiat  plate  endwall,  at  zero  incidence. 
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THE  PRESENT  APPROACH 
Background 

Because  of  the  flow  separation  and  horseshoe  vortex  flow  structure,  analytical 
approaches  based  on  simplifications  such  as  boundary  layer  theory,  rotational  inviscid 
(Euler)  equations,  and  forward  marching  solution  procedures  do  not  seem  adequate  for 
this  problem.  Consequently,  the  approach  being  taken  is  the  numerical  solution  of  the 
compressible  Reynolds-averaged  Navier-Stokcs  equations  using  an  efficient  and  non¬ 
iterative  time-dependent  linearized  block  implicit  (I.B1)  scheme.  A  "zone  embedding" 
approach  is  used  in  which  the  flow  is  computed  only  in  a  subregion  of  the  overall  flow 
field  in  the  immediate  vicinity  of  the  leading  edge.  This  approach  reduces  the  dif¬ 
ficulty  of  constructing  adequate  body-fitted  coordinate  systems  and  is  very  economical., 
although  it  can  he  difficult  to  specify  boundary  conditions  which  adequately  match  the 
embedded  flow  region  to  the  surrounding  flow  region.  A  method  for  applying  inflow/ 
outflow  boundary  conditions  on  curved  coordinate  surfaces  was  developed  in  the  previous 
study  [A]  of  Laminar  horseshoe  vortex  flow  and  is  also  used  here.  These  boundary 
conditions  are  derived  from  an  assumed  flow  structure  and  physical  approximations,  and 
allow  inflow  of  an  inviscid  free  stream  and  developing  boundary  layer,  and  outflow  of 
streamwise  vorticity.  The  boundary  conditions  also  transmit  pressure  waves  through 
the  inflow  boundary  during  the  transient  process,  and  this  avoids  the  instability 
or  slow  convergence  often  attributed  to  t lie  reflection  of  pressure  waves  at  an  inflow 
boundary  with  fixed  velocity  and  pressure.  The  equations  are  solved  in  the  low  Mach 
number  regime  (M  =  0.05)  for  which  they  approximate  the  flow  of  a  liquid.  The  methods 
used  here  have  been  applied  in  the  transonic  and  supersonic  flow  regimes  on  other 
problems . 

Governing  Equations 

The  computation  of  flow  past  swept  leading,  edge  geometries  requires  the  use  of 
nonorthogonal  coordinates.  This  is  accomplished  here  by  .1  general  coordinate  trans¬ 
formation  which  transforms  the  governing  equations  from  a  reference  Cartesian  coordinate 
system  to  a  general  nonorthogonal  coordinate  system  which  is  fitted  to  the  geometry  of 
interest.  The  Cartesian  velcoity  components  u^  are  retained  as  dependent  variables  in 
the  transformed  system  of  equations.  This  same  technique  of  transforming  the  Cartesian 
system  to  general  nonorthogonal  coordinates  is  also  used  by  Pulliam  and  Steger  [b]  and 
Shamroth,  Ctbeling  and  McDonald  (7). 
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The  transformation  T  from  Cartesian  coordinates  to  computational  coordinates 


y-1  is  given  by 


T  =  yJ(x.) 


i,.i  =  1,2,3 


Spatial  derivatives  are  transformed  according  to 

sr-Hirr  <2) 

1  <,y- 

where  unless  otherwise  stated  the  summation  convention  is  used  for  repeated  indices, 
and  y-*  .  =  3yl/3x..  The  coordinate  system  is  defined  by  specifying  the  Cartesian 

» 1  1  i 

coordinates  of  each  computational  grid  point.  The  partial  derivatives  3x^/3y  of  the 
inverse  transformation  T  ^  =  x^(y^)  are  then  computed  using  three-point  second-order 
difference  formulas  with  uniform  spacing  of  the  computational  coordinates  y-* .  The 
transformation  derivatives  3y"*/3x^  are  then  computed  from  3x^/3yJ  using  standard 
procedures  for  computing  derivatives  of  inverse  functions  (cf.  Kaplan  [5]). 

Letting  J  denote  the  Jacobian  determinant  of  the  inverse  transformation 

3 (x. , x^ ,  x~) 

- ± — f.  (3) 

j  =  l  2  3,  v  ; 

3(y  ,y  ,y  ) 

and  making  use  of  the  identity 


the  governing  Navier-Stokes  equations  are  written  in  the  following  nonditnensional 
form:  The  continuity  equation  is 


3(PJ) 

3t 


—  J  yJ  .  pu .  =  0 
J  .  i  i 


where  p  is  density  and  t  is  time. 


The  k-th  component  of  the  momentum  equation  is  given  by 


3(PukJ)  3  •. 

3t  +  J  y',i  (pUiUk  +  6ik  P  "  rik}  =  ° 

where  p  is  pressure  and  the  shear  stress  is  Riven  by 


Tik  =  — 
Re 


(3u. 

m  i  ,  m 

^  a7  y-‘ 


1  3u 

2  .  m  _ l_ 

3  ik  y , i  m 

3y 


Here,  y  is  a  nondimens ional  isotropic  effective  viscosity  coefficient  which  includes 
e 

both  laminar  and  turbulent  contributions.  Re  is  a  reference  Reynolds  number,  and 


b 


is  the  Kroneeker  delta  function. 

The  equation  of  state  for  a  perfect  gas  is  given  hy 

S>  =  pT/yM“  (8) 

where  T  is  temperature,  M  is  a  reference  Mach  number,  and  a  is  the  specific  heat  ratio 

The  total  enthalpy  E  is  defined  by 

E  =  T/(v-l)M“’  +  q  ’  ’  /  -  (9) 

where  q^  =  6*^  u^u  anc^  el  is  t  r  ibn  L  ion  is  governed  by  the  energy  equation. 

Although  it  is  not  necessary,  it  is  both  convenient  and  computationally  worthwhile 
for  the  present  problem  to  assume  thaL  E  is  a  constant,  E  ,  and  to  omit  solution  of 
the  energy  equation.  This  results  in  negligible  error  for  flow  ;il  low  Mach  number 
with  no  heat  addition.  Equations  (8)  and  (b)  can  then  be  combined  to  produce  an 
adiabatic  equation  of  state 

r  =  p(eo  -  q“/2)  (y-D/v  (10) 

which  is  used  to  eliminate  pressure  as  a  dependent  variable  in  Eq .  ((>)  . 

Coordinate  System 

Perspective  views  of  the  geometry,  coordinate  system  and  computational  grid  for 
each  of  the  two  leading-edge  configurations  considered  here  are  shown  in  Fig.  1.  The 
coordinate  system  fits  all  solid  surfaces  within  the  computational  domain  but  is  not 
aligned  with  the  direction  of  the  free  stream  flow.  The  coordinate  system  is  of  suf¬ 
ficient  generality  to  treat  leading  edges  of  arbitrary  (smooth)  cross-sect ional  shape 

and  with  a  swept  leading  edge  whose  angle  may  vary  in  the  spanwise  direction.  In  each 

3  3 

plane  (y  =  constant)  parallel  to  the  endwall,  the  coordinates  consist  of  radial  v 

lines  normal  to  the  curve  where  this  plane  intersects  the  airfoil  surface,  and 

circumferential  y“  lines  parallel  to  this  curve.  These  coordinate  planes  are  "sheared 

by  adding  an  arbitrary  displacement  D(y^)  to  the  Cartesian  x^  coordinates  for  the 

endwall  ;*lane.  Analytical  coordinate  transformations  due  to  Roberts  [81  were  used 

for  each  coordinate  direction  to  redistribute  grid  points  for  adequate  resolution  of 

the  viscous  sublayer  regions  on  both  the  airfoil  and  endwall  surfaces  and  of  the 

geometric  curvature  near  the  leading  edge. 

Turbulence  Model 

The  turbulence  model  used  falls  into  the  category  of  one-equation  turbulence 
models  discussed  by  Launder  and  Spalding  [9],  and  parallels  the  method  given  by 
Shamroth  and  Gibeling  [10].  This  model  requires  solution  of  a  single  partial  dif- 


5 


ferential  equation  governing  turbulence  kinetic  energy  k,  in  eonjuctiou  with  a 

specified  length  scale  i.  A  turbulent  viscosity  p  is  obtained  from  the  Prandtl- 

1/2  t 

Kolmogorov  constitutive  relationship  a  !  k  .  The  turbulent  viscosity  p  is 

assumed  to  be  isotropic,  and  the  stress  tensor  in  the  ensemble-averaged  equation  is 
determined  by  adding  the  turbulent  viscosity  to  the  molecular  viscosity  p  to  obtain 
the  total  effective  viscosity  pg  =  p  +  p  . 

Specification  of  the  length  scale  is  the  remaining  task  which  is  necessary  to 
treat  the  three-dimensional  horseshoe  vortex  flow.  This  flow  has  turbulent  shear 
layers  both  on  the  endwall  (fuselage)  beginning  upstream  of  the  leading  edge  and  also 
on  the  airfoil  or  strut  surface  (and  corner  region)  beginning  at  the  leading  edge. 

Given  an  estimate  of  the  mixing  length  at  the  edge  of  each  of  these  shear  layers  and 
its  streamwise  growth  rate,  a  distribution  of  mixing  length  within  the  shear  layers  can 
be  obtained  from  semi-empirical  relationships  widely  used  in  two-  and  three-dimensional 
boundary  layer  calculations.  To  estimate  the  growth  rate  of  the  two  shear  layers,  a 
boundary  layer  momentum  integral  procedure  is  used  for  each  type  of  shear  layer.  Wind 
tunnel  blockage  effects  are  included  in  the  endwall  shear  layer  calculation,  and 
turbulent  transition  is  presumed  to  occur  at  the  leading  edge  in  the  airfoil  shear 
layer  calculation.  The  final  turbulence  model  provides  for  resolution  of  the  viscous 
sublayer  region  near  walls,  and  permits  calculation  of  the  separated  horseshoe  vortex 
regions  near  the  leading  edge  and  in  the  corner. 

The  equation  governing  turbulence  kinetic,  energy  is  given  by  Launder  and 
Spalding  [11]  for  Cartesian  coordinates.  Applying  the  general  coordinate  transforma¬ 
tion  (1),  this  equation  can  be  written  as 


where 


9 (pkJ) 


at 


3y 


J  y  . 
J  »i 


Re 


,  .  t.  ,  m 

(;,  +  -)  v 

k 


D.  . 
ij 


m 
v  . 
,  i 


(ID 


(12) 


Here,  e  is  the  turbulence  dissipation  rate  and  is  the  turbulent  Prandtl  number 
(taken  here  as  1.0). 


b 


For  high  Reynolds  number  flow,  c  is  approx incite] v  0.09.  For  low  Reynolds  number 

and  in  Che  viscous  sublayer,  e  is  given  a  prescribed  dependence  on  turbulence  Reynolds 

number  R^  in  the  manner  suggested  by  McDonald  and  Fish  [12]  for  transitional  boundary 

layer  flows.  In  [12],  the  turbulence  Reynolds  number  R^  was  appropriately  defined  as  an 

integral  average  across  the  boundary  layer.  Here  and  following  [10],  it  is  convenient 

to  define  R  as  the  local  ratio  of  turbulent  to  laminar  visoositv  R  =  p./p.  The 
T  T  t 

structural  coefficient  c  is  given  in  [12]  us  c  -  i(a.)  ,  where 

M  k*  1 

a.  =  a  f(R  )  /  [1  +  b . 6b  a  (f(R  )-!)]  (15) 

lot  o  r 


with  a  cubic  polynomial  curve  fit  for  values  of  R  between  1  and  40. 

Tt  remains  to  specify  a  length  scale  distribution  appropriate  for  the  problem 
under  consideration.  The  horseshoe-vortex/cornor  I  low  of  interest  here  lias  moderately 
thin  shear  layers  on  the  endwall  and  strut  surf  ice,  and  the  length  scale  distribution 
is  thus  adapted  from  previous  turbulence  models  for  turbulent  boundary  layers  and 
taken  to  be  the  conventionally  defined  mixing  length  of  Prandtl.  The  distribution  of 
mixing  length  given  by  McDonald  and  Fish  [12]  has  proven  effective  for  a  wide  range 
of  two-dimensional  turbulent  boundary  layers  and  is  easily  adapted  for  present  use. 
This  distribution  is  given  by 

i  =  tanh  (ud/f  )  (17) 

i" 


where  £  is  mixing  length,  P.  is  an  outer-region  value  of  mixing  length,  d  is  distance 
from  the  wall,  k  is  the  von  Karman  constant  (taken  here  as  0.41),  and  is  a 

sublayer  damping  function  given  by 

r'/J  [  (d+-2T) /  8 ]  (18) 

+  +  1/2 

Here,  P  is  the  normal  probability  function  and  d  is  defined  by  d  =  d(t/c)  /v, 
where  t  is  shear  stress  and  v  is  kinematic  viscosity.  For  equilibrium  boundary  layers. 
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Z  is  observed  to  have  a  constant  value  of  about  0.09  6,  where  6  is  the  local  boundary 
layer  thickness. 

The  length  scale  distribution  of  Eq.  (17)  is  adapted  for  present  use  by  taking 
d  as  distance  to  the  nearest  wall  and  by  assigning  f  a  distribution  based  on  two- 
dimensional  momentum  integral  estimates  of  the  boundary  layer  growth  rates  on  the 
endwall  and  on  the  strut.  The  computed  estimates  of  boundary  layer  growth  were 
obtained  from  a  simple  integral  prediction  scheme  rather  than  an  attempt  to  scan  the 
intermediate  transient  solutions  of  the  Navier-Stokes  equations  to  determine  some  ill- 
defined  boundary  layer  thickness.  Assuming  a  1/7  power  lav  velocity  profile  and  a 
skin  friction  law  c^/2  =  0.0225  (v/u  the  momentum  integral  equation  can  be 

written  as 

, ,  .  du 

-  V2  <19> 

e 

where  for  the  1/7  power  profile:  h^  =  7/72  and  =  23/72.  The  momentum  integral 

equation  is  applied  to  both  strut  and  endwall  boundary  layers.  The  free  stream 

velocity  ug  is  imposed  from  a  potential  flow.  The  flow  cases  computed  here  are 

enclosed  in  a  wind  tunnel,  and  since  the  endwall  boundary  layer  causes  blockage 

effects  in  this  instance,  the  free  stream  velocity  is  adjusted  for  this  effect.  If 

distances  are  normalized  by  the  tunnel  half-height,  then  the  free  stream  velocity  adjusted 

for  blockage  is  given  by  u  =  1 / (1—6  ),  where  for  the  1/7  power  velocity  profile,  the 

*  e  * 

displacement  thickness  <5  is  related  to  5  by  6  =  8  6  .  This  relationship  for  ug  is 
inserted  into  Eq.  (19)  prior  to  solution.  Eq.  (19)  is  solved  using  a  second-order 
linearized  finite  difference  scheme  described  in  [13] . 

Two  boundary  layer  thickness  distributions,  6.  (x. )  for  the  endwall  boundary 
2 

layer  and  6 ^  (y  )  for  the  strut  surface,  are  obtained  by  solution  of  Eq.  (19).  An 
outer  or  maximum  mixing  length  scale  associated  with  each  boundary  layer  is  determined 
from  the  formula  Z  =  0.09  6.  In  the  core  region  outside  the  boundary  layer,  the 
length  scale  is  exponentially  damped  to  a  negligible  value  over  a  distance  of  about 
two  boundary  layer  thicknesses.  In  the  overlap  corner  region,  the  outer  or  maximum 
length  scale  Z  is  taken  as  the  greater  of  the  two  values  associated  with  the  strut  and 
endwall  boundary  layers. 

Physical  Boundary  and  Initial  Conditions 

Since  the  computational  domain  is  chosen  to  be  a  region  in  the  immediate  vicinity 
of  the  leading-edge/corner  flow  geometry  (cf.  Fig.  1)  embedded  within  a  larger  overall 
flow  system,  inflow  and  outflow  boundary  conditions  which  adequately  model  the  inter¬ 
face  between  the  computed  flow  and  the  remainder  of  the  flow  system  are  required. 
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The  inflow/outflow  conditions  used  are  derived  from  an  assumed  flow  structure  and  are 
chosen  to  provide  inflow  with  prescribed  stagnation  pressure  (and  stagnation  enthalpy) 
in  an  inviscid  core  region  and  with  a  given  axial  velocity  profile  shape  in  the  endwall 
boundary  layer,  and  to  provide  outflow  with  a  prescribed  distribution  of  static  pres¬ 
sure  in  the  cross  section.  These  boundary  conditions  are  compatible  both  with  an 
inviscid  characteristics  analysis  and  with  the  physical  process  by  which  most  flows 
are  established.  These  boundary  conditions  allow  both  velocity  and  static  pressure  to 
vary  with  time  at  the  inflow  boundary,  and  consequently,  pressure  waves  are  transmitted 
upstream  through  the  inflow  boundary  during  the  transient  flow  process  and  are  not 
reflected  back  into  the  computational  domain.  The  reflection  of  pressure  waves  at  an 
inflow  boundary  where  velocity  and  pressure  are  fixed  in  time  has  often  been  cited  as 
a  cause  of  either  instability  or  slow  convergence  in  other  investigations.  These 
boundary  conditions  are  discussed  in  more  detail  in  [14], 

The  initial  and  boundary  conditions  are  devised  from  a  rough  approximation  of  the 
two-dimensional  potential  flow  velocity  for  the  strut  cross  section,  and  from  the 

momentum-integral  estimates  for  blockage  B(x.)  and  for  the  boundary  l.ayer  thickness 

2  1 

distributions  S^(x^)  and  6 2 (y  )  on  the  endwall  flat  plate  and  strut  surfaces,  respectively. 

Finally,  boundary  layer  velocity  profile  shapes  f^  (y/6^)  and  f^  (y/^),  0  -  f^, 

f2  -  1  are  defined,  where  y  is  a  parameter  indicative  of  distance  from  a  wall.  These 

velocity  profile  shapes  were  taken  from  the  analytical  fit  of  Musker  [15]  to  the  Coles- 

* 

type  of  profile  which  matches  6  and  c^.  from  the  moment urn- integral  calculations.  The 
initial  conditions  at  time  t  =  o  are  defined  by 

U  =  UT  B(xx)  f1  (y/6x)  f2(y/«2)  (20) 

c  =  1  -  B2  UT  •  UT  (21) 

P  I  I 

where  U  is  the  velocity  vector  and  c^  is  the  pressure  coefficient  referred  to  the 
reference  conditions.  The  details  of  this  procedure  are  not  critical  and  are  omitted, 
since  except  for  the  shape  and  thickness  of  the  inlet  boundary  layer  profile,  these 
results  serve  only  as  a  convenient  method  for  selecting  initial  conditions.  A  reasonably 
accurate  estimate  for  the  pressure  drop  which  will  produce  the  desired  flow  rate  must 
be  made  using  any  convenient  source.  For  the  present  problem,  the  approximate  poten¬ 
tial  flow  corrected  for  blockage  as  in  Eq.  (21)  is  adequate.  It  is  noted  that 
although  these  initial  conditions  do  take  into  account  several  relevant  features 
of  the  flow,  the  important  effects  of  flow  separation,  and  horseshoe/corner  vortex 
formation  are  completely  neglected.  The  initial  flow  field  is  thus  a  simple  but 
relatively  crude  approximation  to  the  final  flow  field. 


a  "two-layer"  boundary  condition  is  employed 


r 

p 


At  the  inflow  boundary  y^  = 

as  in  [4],  such  that  stagnation  pressure  p  is  fixed  at  the  free  stream  reference 

3  0 

value  in  the  core  flow  region  (y  >  6,)  and  the  axial  velocity  profile  shape  u./u  = 

_  -*•  2  i  e 

f^(y/6^)  is  fixed  within  the  boundary  layer  region  (y  -  6^).  Here,  ug  is  the  local 

edge  velocity  which  varies  with  time  and  is  adjusted  after  each  time  step  to  the  value 

consistent  with  p  and  the  local  edge  static  pressure,  which  is  determined  as  part  of 

°  2  2  2  2 
the  solution.  The  remaining  inflow  conditions  are  U2  =  3  u^/Sn  =  3  c^/3n  = 

where  n  denotes  the  normal  computational  coordinate,  y^.  For  outflow  conditions,  the 

static  pressure  is  imposed,  and  second  derivatives  of  each  velocity  component  are  set 

to  zero. 

At  no-slip  surfaces,  each  velocity  component  u^  is  set  to  zero,  and  the  remaining 
condition  applied  at  these  surfaces  is  that  the  derivative  of  pressure  in  the  direction 
normal  to  the  surface  is  zero.  This  condition  approximates  the  normal  momentum  equa¬ 
tion  to  order  Re  ^  for  viscous  flow  at  a  no-slip  surface.  For  the  swept  leading  edge 
configuration,  the  surface  normal  does  not  always  lie  on  a  coordinate  line,  and  the 
normal  derivative  boundary  condition  must  be  derived  from  the  coordinate  transforma¬ 
tion  data.  Letting  e.  =  3r/3y^  denote  the  basis  vectors  for  y^ ,  where  r  is  the 

I  _  12 

position  vector,  then  the  unit  vector  n  normal  to  the  y  -  y  coordinate  surface  is 
given  by  n  =  e^  x  e2/lei  x  e2 I •  If  n*  are  the  Cartesian  components  of  n,  then  the 


normal  derivative  of  <j>  is  given  by 


-  i  j  3d> 

n  •  V<J>  =  n  y  .  — . 

,X  3yJ 


(22) 


The  final  boundary  to  be  considered  is  the  plane  parallel  to  the  endwall  and 
in  the  free  stream.  This  boundary  is  assumed  to  be  a  plane  of  symmetry,  so  that  the 
flow  represented  is  that  past  the  strut  mounted  between  parallel  flat  plates  with 
spacing  2H.  This  assumption  is  both  convenient  and  corresponds  to  the  conditions  in  the 
wind  tunnel  experiment  for  the  unswept  leading  edge  flow,  for  which  experimental 
measurements  are  available. 


SOLUTION  PROCEDURE 
Background 


The  solution  procedure  employs  a  consistently-split  linearized  block  implicit 
(LBI)  algorithm  which  has  been  discussed  in  detail  by  the  authors  in  [13,  16].  There 
are  two  important  elements  of  this  method: 

(1)  the  use  of  a  noniterative  formal  time  linearization  to 
produce  a  fully-coupled  linear  multidimensional  scheme 
which  is  written  in  "block  implicit"  form;  and 

(2)  solution  of  this  linearized  coupled  scheme  using  a  consistent 
"splitting"  (ADI  scheme)  patterned  after  the  Douglas-Gunn  [17] 

(1964)  treatment  of  scalar  ADI  schemes. 

The  method  is  thus  referred  to  as  a  split  linearized  block  implicit  (LBI)  scheme. 

The  method  has  several  attributes: 

(1)  the  noniterative  linearization  is  efficient; 

(2)  the  fully-coupled  linearized  algorithm  eliminates  instabilities 
and/or  extremely  slow  convergence  rates  often  attributed  to 
methods  which  employ  ad  hoc  decoupling  and  linearization  assumptions 
to  identify  nonlinear  coefficients  which  are  then  treated  by  lag  and 
update  techniques; 

(3)  the  splitting  or  ADI  technique  produces  an  efficient  algorithm 
which  is  stable  for  large  time  steps  and  also  provides  a  means 
for  convergence  acceleration  for  further  efficiency  in  computing 
steady  solutions; 

(4)  intermediate  steps  of  the  splitting  are  consistent  with  the  governing 
equations,  and  this  means  that  the  "physical"  boundary  conditions  can 
be  used  for  the  intermediate  solutions.  Other  splittings  which  are 
inconsistent  can  have  severe  difficulties  in  satisfying  physical 
boundary  conditions  [13]. 

(5)  the  convergence  rate  and  overall  efficiency  of  the  algorithm  are  much 
less  sensitive  to  mesh  refinement  and  redistribution  than  algorithms 
based  on  explicit  schemes  or  which  employ  ad  hoc  decoupling  and  lineari¬ 
zation  assumptions.  This  is  important  for  accuracy  and  for  computing 
turbulent  flows  with  viscous  sublayer  resolution;  and 
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(6)  the  method  is  general  and  is  specifically  designed  for  the  complex 
systems  of  equations  which  govern  viscous  flow  in  complicated 
geometries. 

This  same  algorithm  was  later  considered  by  Beam  and  Warming  [18],  but  the  ADI 
splitting  was  derived  by  approximate  factorization  instead  of  the  Douglas-Gunn 
procedure.  They  refer  to  the  algorithm  as  a  "delta  form"  approximate  factorization 
scheme.  This  scheme  replaced  an  earlier  non-delta  form  scheme  [19],  which  has 
inconsistent  intermediate  steps. 


Spatial  Differencing  and  Artificial  Dissipation 


The  spatial  differencing  procedures  used  are  a  straightforward  adaptation  of 
those  used  in  [16]  and  elsewhere.  Three-point  central  difference  formulas  are  used 
for  spatial  derivatives,  including  the  first-derivative  convection  and  pressure 
gradient  terms.  This  has  an  advantage  over  one-sided  formulas  in  flow  calculations 
subject  to  "two-point"  boundary  conditions  (virtually  all  viscous  or  subsonic  flows), 
in  that  all  boundary  conditions  enter  the  algorithm  implicitly.  In  practical  flow 
calculations,  artificial  dissipation  is  usually  needed  and  is  added  to  control  high- 
frequency  numerical  oscillations  which  otherwise  occur  with  the  central-difference 
formula. 

In  the  present  investigation,  artificial  (anisotropic)  dissipation  terms  of 
the  form 


a2 

T  i  2  3  uk 

.  .  d  .  J  (y 1  . )  Z  - ^ r 

X’J  J  ;)rV 


(23) 


are  added  to  the  right-hand  side  of  each  (k-th)  component  of  the  momentum  equation  (6), 
where  for  each  coordinate  direction  y1,  the  artificial  diffusivity  d.  is  positive  and 


is  chosen  as  the  larger  of  zero  and  the  local  quantity  p  (a  Re^  j 


the  local  cell  Reynolds  number  Re 


Ayl 


1)/Re.  Here, 
for  the  j-th  direction  is  defined  by 


Re  .  =  Re  | T  y' .  pu.|  AyJ/p  (y1.)2  (24) 

AyJ  i  ’L  1  e  i  -.i 

with  no  summation  on  j.  This  treatment  lowers  the  formal  accuracy  to  0  (Ax),  but  the 
functional  form  is  such  that  accuracy  in  representing  physical  shear  stresses  in  thin 
shear  layers  with  small  normal  velocity  is  not  seriously  degraded.  This  latter 
property  follows  from  the  anisotropic  form  of  the  dissipation  and  the  combination  of 
both  small  normal  velocity  and  small  grid  spacing  in  thin  shear  layers.  A  value  of 
0.5  was  used  for  a  in  the  present  calculations.  Values  lower  than  0.5  have  been  used 
to  good  effect  in  two  space  dimensions  [20,  21],  but  it  has  not  yet  been  possible  to 
investigate  the  role  of  smaller  o  values  in  three  space  dimensions. 


Split  LBI  Algorithm 


Linearization  and  Time  Differencing 

The  system  of  governing  equations  to  be  solved  consists  of  continuity  (5),  three 
components  of  momentum  (6)  and  the  turbulence  kinetic  energy  equation  (11).  Ancillary 
relations  (10)  and  ( 1 A )  are  substituted  in  these  equations  to  eliminate  p  and  e  as 
dependent  variables.  This  produces  a  system  of  five  equations  in  five  dependent 
variables:  p,  u^,  u3>  and  k.  Using  notation  similar  to  that  in  [16],  at  a 

single  grid  point  this  system  of  equations  can  be  written  in  the  following  form: 


3  H(<(>) /3t  =  D«>)  +  S(4>) 


(25) 


*1 


where  (J>  is  the  column-vector  of  dependent  variables,  H  and  S  are  column-vector 
algebraic  functions  of  <J>,  and  D  is  a  column  vector  whose  elements  are  the  spatial 
differential  operators  which  generate  all  spatial  derivatives  appearing  in  the  govern¬ 
ing  equation  associated  with  that  element.  In  the  present  case,  the  only  non-zero 
i  element  in  S  is  generated  by  the  last  term  in  the  turbulence  kinetic  energy  equation 

(11). 

The  solution  procedure  is  based  on  the  following  two-level  implicit  time- 
difference  approximations  of  (25): 

(Hn+1  -  Hn)/At  =  6(Dn+1  +  Sn+1)  +  (1-6)  (Dn  +  Sn)  (26) 

where,  for  example,  Hn+^  denotes  H((j>n+'*')  and  At  =  tn+^  -tn.  The  parameter  B 

(0.5  -3-1)  permits  a  variable  time-centering  of  the  scheme,  with  a  truncation 
2 

error  of  order  [At  ,  (6  -  1/2)  At]. 

A  local  time  linearization  (Taylor  expansion  about  4>n )  of  requisite  formal 
accuracy  is  introduced,  and  this  serves  to  define  a  linear  differential  operator  L 
(cf.  [16])  such  that 

Dn+1  =  Dn  +  Ln  (<J>n+1  -  <j>n)  +  0  (At2)  (27) 

Similarly, 

Hn+1  =  Hn  +  (3H/3<f>)n  (<i>n+1  -  $n)  +  0  (At2)  (28) 

Sn+1  =  Sn  +  (3S/3<f>)n  (<f>n+1  -  ]>n)  +  0  (At2)  (29) 

Eqs.  (27-29)  are  inserted  into  Eq.  (26)  to  obtain  the  following  system  which  is 
linear  in  <f>n+^ : 

(A  -  BAt  Ln)  (<{,n+1  -  (f.n)  =  At  (Dn  +  Sn)  (30) 

and  which  is  termed  a  linearized  block  implicit  (LBI)  scheme. 
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Here,  A  denotes  a  square  matrix  defined  by 


A  =  (<)H/0cf>)n  -  3At  (3S/3<J>) 


(31) 


Eq.  (30)  has  0  (At)  accuracy  unless  H  =  <j>,  in  which  case  the  accuracy  is  the  same  as 
Eq.  (26). 

Special  Treatment  of  Diffusive  Terms 

The  time  differencing  of  diffusive  terms  is  modified  to  accomodate  cross¬ 
derivative  terms  and  also  turbulent  viscosity  and  artificial  dissipation  coefficients 
which  depend  on  the  solution  variables.  Although  formal  linearization  of  the  convection 
and  pressure  gradient  terms  and  the  resulting  implicit  coupling  of  variables  is  critical 
to  the  stability  and  rapid  convergence  of  the  algorithm,  this  does  not  appear  to  be 
important  for  the  turbulent  viscosity  and  artificial  dissipation  coefficients.  Since 
the  relationship  between  and  d^  and  the  mean  flow  variables  is  not  conveniently 
linearized,  these  diffusive  coefficients  are  evaluated  explicitly  at  tn  during  each 
time  step.  Notationally ,  this  is  equivalent  to  neglecting  terms  proportional  to 
3pe/3(J>  or  9dj/3<|>  in  Ln,  which  are  formally  present  in  the  Taylor  expansion  (27),  but 
retaining  all  terms  proportional  to  or  d^  in  both  Ln  and  Dn. 

It  has  been  found  through  extensive  experience  that  this  has  little  if  any 
effect  on  the  performance  of  the  algorithm.  This  treatment  also  has  the  added 
benefit  that  the  turbulence  model  equations  (in  this  instance  the  turbulence  kinetic 
energy  equation)  can  be  decoupled  from  the  system  of  mean  flow  equations  by  an  appro¬ 
priate  matrix  partitioning  (cf.  [13])  and  solved  separately  in  each  step  of  the  ADI 
solution  procedure.  This  reduces  the  block  size  of  the  block  tridiagonal  systems 
which  must  be  solved  in  each  step  and  thus  reduces  the  computational  labor. 

In  addition,  the  viscous  terms  in  the  present  formulation  include  a  number  of 
spatial  cross-derivative  terms.  Although  it  is  possible  to  treat  cross-derivative 
terms  implicitly  within  the  ADI  treatment  which  follows,  it  is  not  at  all  convenient 
to  do  so,  and  consequently,  all  cross-derivative  terms  are  evaluated  explicitly  at  tn. 
For  a  scalar  model  equation  representing  combined  convection  and  diffusion,  it  has 
been  shown  by  Beam  and  Warming  [22]  that  the  explicit  treatment  of  cross-derivative 
terms  does  not  degrade  the  unconditional  stability  of  the  present  algorithm.  To 
preserve  notational  simplicity,  it  is  understood  that  all  cross-derivative  terms 
appearing  in  Ln  are  neglected  but  are  retained  in  Dn.  It  is  important  to  note  that 
neglecting  terms  in  Ln  has  no  effect  on  steady  solutions  of  Eq.  (30),  since 
$n+^-<j>n  =  0  and  thus  Eq.  (30)  reduces  to  the  steady  form  of  the  equations:  D°  +  Sn  =  0. 
Aside  from  stability  considerations,  the  only  effect  of  neglecting  terms  in  Ln  is  to 
introduce  an  0  (At)  truncation  error. 


Consistent  Splitting  of  the  LBI  Scheme 


To  obtain  an  efficient  algorithm,  the  linearized  system  (30)  is  split  using 
ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional  operator  L  is 
rewritten  as  the  sum  of  three  "one-dimensional"  sub-operators  (i  =  1,  2,  3)  each  of 
which  contains  all  terms  having  derivatives  with  respect  to  the  i-th  spatial 
coordinate.  The  split  form  of  Eq.  (30)  can  be  derived  either  as  in  [13,  16]  by 
following  the  procedure  described  by  Douglas  and  Gunn  [17]  in  their  generalization  and 
unification  of  scalar  ADI  schemes,  or  using  approximate  factorization  as  in  [18], 

For  the  present  system  of  equations,  the  split  algorithm  is  given  by 


(A  - 

BAtL^) 

* 

(♦  - 

*n)  = 

At 

(Dn  + 

sn) 
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where  <|>  and  (f>  are  consistent  intermediate  solutions  [13,  16].  If  spatial  deriva¬ 
tives  appearing  in  and  D  are  replaced  by  three-point  difference  formulas,  as 
indicated  previously,  then  each  step  in  Eqs.  (32a-c)  can  be  solved  by  a  block- 
tridiagonal  elimination. 

Combining  Eqs.  (32a-c)  gives 

(A  -  BAtL")  A_l  (A  -  BAtL^)  A_1  (A  -  BAtL^)  (4>n+1  -  <j>n) 

n  (33) 

=  At  (D  +  S  ) 

2 

which  approximates  the  unsplit  scheme  (30)  to  0  (At  ).  Since  the  intermediate  steps 
are  also  consistent  approximations  for  Eq.  (30),  physical  boundary  conditions  can  be 

'It  Idt 

used  for  <|>  and  <j>  [13,  16].  Finally,  since  the  are  homogeneous  operators, 

it  follows  from  Eqs.  (32a-c)  that  steady  solutions  have  the  property  that 

n+1  *  **  n 

d>  =  and  satisfy 

Dn  +  Sn  =  0  (34) 


The  steady  solution  thus  depends  only  on  the  spatial  difference  approximation  used 
for  (34)  and  does  not  depend  on  the  solution  algorithm  itself. 


COMPUTED  RESULTS 


Solutions  for  three-dimensional  turbulent  flow  past  both  unswept  and  45-degree 
swept  elliptical  leading  edges  affixed  to  a  flat  plate  endwal 1  are  presented  here. 

To  assess  the  degree  of  accuracy  obtainable  in  three-dimensional  flow  calculations,  it 
is  obviously  valuable  to  have  experimental  measurements  for  comparison.  The  recent 
measurements  by  Mehta,  Shabaka  and  Bradshaw  [1]  for  turbulent  corner  flow  downstream 
of  an  unswept  elliptical  leading  edge  at  zero  incidence  appear  to  be  the  only  measure¬ 
ments  available  which  include  secondary  velocity  components,  and  this  flow  was  selected 
for  computation. 

Solutions  were  computed  for  turbulent  horseshoe  vortex  flow  past  both  unswept 

and  45-degree  swept  leading  edges  having  the  configurations  shown  in  Fig.  1,  each 

having  a  Reynolds  number  Re  =  130,000  and  Mach  number  M  =  0.05.  The  grids  used  were 

12  3 

18  x  28  x  18  and  18  x  24  x  21  (y  ,  y  ,  y  )  for  the  unswept  and  swept  cases,  respectively. 
The  Reynolds  number  is  based  on  the  tunnel  half-height  of  2.5  inches.  The  geometry 
consists  of  an  elliptical  leading  edge  with  a  half-chord  of  6  inches,  followed  by  a 
slab  with  constant  thickness  of  2  inches.  In  the  configuration  with  swept  leading 
edge,  this  cross  section  is  "sheared”  as  shown  in  Fig.  l.b  to  form  a  sweep  angle  of 
45  degrees  at  the  endwall.  The  angle  of  sweep  varies  smoothly  between  45  degrees 
at  the  endwall  and  zero  degrees  at  the  symmetry  plane  midway  between  the  wind 
tunnel  walls.  This  particular  swept  configuration  was  selected  to  produce  a  flow 
which  is  symmetric  about  this  midplane,  for  ease  in  imposing  boundary  conditions. 

Mehta,  Shabaka  and  Bradshaw  indicate  that  the  boundary  layer  on  the  endwal 1  has  a 
thickness  of  approximately  1.0  inches  at  the  leading  edge;  momentum  integral  cal¬ 
culations  were  used  to  estimate  an  inflow  boundary  layer  thickness  which  would  match 
this  condition. 

Convergence  Rate 

Turbulent  flows,  particularly  in  three  dimensions,  are  especially  difficult 
to  compute  because  of  the  very  high  degree  of  mesh  redistribution  required  to  resolve 
viscous  sublayer  regions.  Resolution  of  the  viscous  sublayer  is  believed  necessary 
to  obtain  an  adequate  representation  of  the  turbulent  shear  layer  behavior.  In 
addition,  because  of  the  leading-edge  geometry  and  the  high  strut  length-to-thickness 
aspect  ratio,  this  mesh  redistribution  is  required  for  all  three  coordinate  directions. 
This  contrasts,  for  example,  with  three-dimensional  flow  in  a  curved  pipe  or  duct, 
where  only  one  or  two  (respectively)  coordinate  directions  require  substantial  mesh 
redistribution.  As  a  consequency,  the  ratio  of  largest  to  smallest  elemental 
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control  volume  for  the  present  calculation  is  approximately  3  x  10  .  The  correspond¬ 
ing  matrix  of  difference  approximations  has  an  extremely  wide  range  of  eigenvalues, 
and  thus  poses  a  computational  problem  for  which  rapid  convergence  is  difficult  to 
achieve . 

The  present  computational  method  employs  variable  time  steps  to  accelerate 
convergence,  as  discussed  in  [23].  Several  procedures  for  time  step  selection  are 
presently  under  evaluation,  but  typically  small  time  steps  are  used  initially  during 
the  elimination  of  impulsive  transients,  followed  by  larger  time  steps  and  the  cyclic 
use  of  a  sequence  of  time  steps.  In  addition,  a  smaller  time  step  was  used  near  the 
leading  edge  than  elsewhere  in  the  flow  field.  Small  time  steps  are  effective  in 
reducing  high  (spatial)  frequency  error  components,  while  large  time  steps  are 
effective  in  reducing  low  frequency  errors. 

An  indication  of  both  the  degree  and  rate  of  convergence  obtained  for  the 
present  calculations  is  shown  in  Fig.  2.  For  comparison,  a  curve  representative  of 
the  convergence  behavior  typically  obtained  for  two-dimensional  turbulent  flow  cases 
is  also  shown.  Although  the  present  convergence  rate  is  noticeably  slower  than  that 
observed  using  this  same  algorithm  in  other  calculations,  it  is  adequate  for  present 
purposes.  The  present  calculations  ere  terminated  after  approximately  100  iterations 
for  reasons  of  economy,  since  it  appeared  that  changes  in  the  solution  at  this  point 
were  of  minor  significance  and  would  not  alter  any  of  the  conclusions  drawn  from 
these  results. 

Results  for  Unswept  Leading  F.dge 

Computed  results  for  the  unswept  leading-edge  case  are  given  in  a  series  of 
plots  in  Figs.  3  -  10.  It  is  difficult  to  display  grapnicallv  the  results  of  three- 
dimensional  flow  calculations,  and  the  present  approach  is  to  consider  contour  and 
velocity  vector  plots  for  selected  two-dimensional  planes  and  projected  surfaces, 
which  are  representative  of  the  computed  flow  structure.  Because  of  the  multiple 
length  scales  in  these  results,  the  solutions  are  plotted  in  some  instances  for  the 
entire  computational  region  and  grid,  and  in  other  instances  only  a  portion  of  the 
computed  solution  is  shown.  Grid  point  indices  are  used  in  these  plots  as  a  means  of 
identifying  the  extent  of  the  region  shown.  The  notation  x,  y,  z  is  used  to  denote 
the  respective  computational  coordinates  y1  shown  in  Fig.  1 . 

In  Fig.  3,  the  coordinates  and  computational  grid  are  shown  for  planes  which 
are  parallel  to  the  endwall.  In  Figs.  4  and  5,  velocity  vector  plots  from  the 
nominally  two-dimensional  flow  region  at  the  free  stream  symmetry  plane  midway  between 
the  tunnel  walls  are  contrasted  with  the  corresponding  results  at  the  plane  of  grid 
points  adjacent  to  the  endwall  and  very  near  the  endwall  surface.  The  vector 
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magnitudes  are  renormalized  for  each  plot  and  thus  indicate  only  flow  direction  and 
relative  magnitude  within  the  plot.  In  this  region  near  the  leading  edge,  the 
computed  flow  structure  has  all  of  the  characteristics  observed  in  numerous  flow 
visualization  experiments  and  in  the  authors'  previous  laminar  calculations  f4J. 

The  flow  near  the  endwall  separates  and  forms  a  horseshoe  vortex  upstream  of  the 
leading  edge  with  peak  reverse  flow  velocity  in  this  instance  about  20  per  cent 
of  the  free  stream  velocity.  The  reversed  flow  within  the  horseshoe  vortex  is 
clearly  evident  in  Fig.  5.  A  secondary  vortex  can  also  be  seen  at  the  leading  edge, 
but  this  is  only  marginally  resolved. 

In  Fig.  6,  the  static  pressure  and  velocity  are  shown  for  the  (vertical)  stagna¬ 
tion  flow  plane  containing  the  leading  edge  and  extending  upstream  parallel  to  the 
freestream  flow.  The  horseshoe  vortex  is  too  small  and  too  close  to  the  leading 
edge  to  be  visible  on  the  scale  of  Fig.  6,  but  is  shown  in  Fig.  7  in  a  detail  of  the 
region  near  the  intersection  of  the  leading  edge  and  endwall.  The  peak  velocity  is 
directed  toward  the  endwall  and  parallel  to  the  leading  edge,  and  is  about  35  per  cent 
of  the  freestream  velocity. 

A  full  scale  view  of  the  axial  velocity  and  secondary  or  crossflow  velocity  is 
shown  in  Fig.  8  for  the  plane  perpendicular  to  the  corner  intersection  and  located 
at  the  end  of  the  elliptical  leading  edge  region  (the  y  =  20  plane,  cf.  Fig.  3). 
Details  of  the  corner  flow  downstream  of  the  leading  edge  are  shown  in  Figs.  9  and 
10  at  each  of  three  axial  locations  (see  Fig.  3  for  locations).  The  growth  of  the 
corner  shear  layer  and  a  rapid  decay  of  the  secondary  flow  strength  can  be  seen  in 
Figs.  9  and  10.  The  last  axial  location  shown  (y  =  26)  corresponds  approximately  to 
axial  station  2  (Fig.  la)  for  which  measurements  of  both  axial  and  transverse 
velocities  were  taken  by  Mehta,  Shabaka  and  Bradshaw  [1],  Their  measurements  are 
shown  in  Fig.  11,  and  in  this  region  (about  25  inches  downstream  of  the  leading  edge) 
the  measurements  show  a  single  corner  vortex  having  the  same  sense  as  that  computed 
near  the  leading  edge  but  with  a  peak  cross  flow  velocity  of  about  5  per  cent  of  the 
free  stream  velocity.  The  computed  cross  flow  does  not  agree  well  with  measurements 
in  this  region,  however,  and  indicates  that  the  corner  vortex  flow  has  very  small 
cross  flow  velocities  (less  than  0.5  per  cent  of  free  stream  velocity),  and  of  op¬ 
posite  sense.  This  poor  agreement  is  presumably  the  result  of  numerical  truncation 
error  and/or  inadequacy  of  the  turbulence  model .  It  is  not  possible  to  make  any 
further  quantitative  assessment  on  the  basis  of  the  available  measurements  alone. 

It  is  unfortunate  that  experimental  measurements  were  not  taken  near  the  leading 
edge,  and  thus  comparisons  cannot  be  made  in  this  critical  region  where  the  vortex 
is  formed. 


Related  difficulties  have  been  encountered  in  obtaining  accurate  three- 
dimensional  flow  predictions  for  flow  in  curved  ducts  [24],  and  in  that  study,  the 
turbulence  model  does  not  appear  to  be  responsible.  It  Is  believed  that  the  present 
disagreement  with  experiment  is  most  likely  the  result  of  numerical  error  including 
the  added  artificial  dissipation.  In  a  recent  invest igat ion  of  two-dimensional  viscous 
transonic  flow  past  airfoils  120],  it  was  found  that  significant  reductions  in  the 
magnitude  (factors  of  10  to  20)  of  the  artificial  dissipation  are  possible  without 
destabilizing  the  calculation  and  with  a  corresponding  increase  in  accuracy.  The 
present  calculation  was  repeated  with  a  tenfold  reduction  in  dissipation  parameter 
but  was  found  to  be  unstable.  Further  work  is  needed  to  resolve  the  present 
discrepancy  between  prediction  and  experiment. 

Results  for  Swept  Leading  Edge 

Computed  results  for  the  45-degree  swept  leading  edge  configuration  in  Fig.  lb 
are  given  in  Figs.  12  -  21.  This  flow  also  separates  and  forms  a  horseshoe  vortex 
upstream  of  the  leading  edge.  The  computed  flow  structure  for  the  swept  geometry 
differs  from  that  for  the  unswept  geometry  in  two  important  respects.  First,  the 
strength  of  the  horseshoe  vortex  at  the  leading  edge  is  significantly  reduced  bv  the 
swept  leading  edge,  as  evidenced  by  a  peak  reverse  flow  velocity  of  approximately 
10  per  cent  of  the  free  stream  velocity  (compared  with  20  per  cent  for  the  unswept 
geometry).  Secondly,  outside  the  horseshoe  vertex  region,  a  cross  flow  parallel  to 
the  swept  leading  edge  and  directed  outward  aw.iv  from  the  endwall  occurs  due  to  the 
swept  geometry.  The  peak  velocity  parallel  to  the  swept  leading  edge  is  predicted 
to  be  approximately  40  per  cent  of  the  free  stream  velocity.  This  cross  flow  dis¬ 
appears  downstream  of  the  swept  leading  edge  and  is  completely  absent  in  the  unswept 
geometry. 

The  grid  distribution  used  for  this  case  is  shown  in  Fig.  12.  Side-view 
projections  of  the  flow  near  the  surface  of  Lite  strut  are  shown  in  Fig.  13  (entire 
grid)  and  in  Fig.  14  (detail  near  leading  edge).  The  large  cross  flow  parallel  to 
the  swept  leading  edge  is  evident  in  each  of  these  figures.  Computed  results  for 
the  stagnation  plane  upstream  of  the  leading  edge  are  shewn  in  Fig.  15  (entire  grid) 
and  Fig.  16  (detail  near  intersection  of  leading  edge  and  endwall).  Computed  dis¬ 
tributions  of  velocity  and  pressure  are  given  in  Fig.  16  for  a  plane  parallel  to  the 
flat  plate  and  located  near  the  edge  of  the  flat  plate  boundary  layer.  Details  of 
the  flow  near  the  leading  edge  are  shown  in  Figs.  IS  and  19  for  the  free  stream 
symmetry  plane  and  for  the  plane  of  grid  points  adjacent  to  the  endwall  surface. 

The  reversed  flow  region  inside  the  horseshoe  vortex  is  clearly  visible  in  Fig.  19. 
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Finally,  projections  of  the  corner  flow  downstream  of  the  l.ading  edge  are  shown 
in  Figs.  20  and  21  at  each  of  three  axial  locations  (see  Fig.  12  for  location). 

To  summarize  these  results,  the  computed  flow  field  near  the  swept  leading  edge 
displays  a  reduced  vortex  strength,  strong  cross  flow  along  the  swept  leading  edge, 
and  a  modified  pressure  distribution  near  the  strut-endwall  intersection.  Although 
these  flow  characteristics  are  qualitatively  what  one  expects  for  this  flow,  no  other 
analytical  or  experimental  results  are  available  for  quantitative  evaluation,  and  at 
this  point  the  results  serve  mainly  as  a  demonstration  of  capability  lor  computing 
this  complex  and  difficut  three-dimensional  turbulent  horseshoe  vortex  flow. 

Concluding  Remarks 

Three-dimensional  turbulent  horseshoe  vortex  flow  past  both  swept  and  unswept 
leading-edge/endwall  configurations  has  been  studied  by  numerical  solution  of  the 
Reynolds-averaged  Navier-Stokes  equations.  It  is  believed  that  the  physical  processes 
involved  require  the  viscous  sublayer  to  be  resolved,  and  the  computational  approach 
provides  for  this  resolution  and  has  been  used  with  general  nonorthogonal  coordinates. 

Predictions  of  the  horseshoe  vortex  formation  near  the  leading  edge  are  in 
qualitative  agreement  with  flow  visualization  studies  of  similar  flows.  No  flow 
measurements  are  available  for  the  region  near  the  leading  edge.  Predictions  of  Ihe 
relatively  weak  secondary  flows  in  the  corner  downstream  of  the  unswept  leading  edge 
do  not  agree  with  available  measurements  well  downstream  of  the  leading  edge.  This 
disagreement  is  believed  to  be  mainly  the  result  of  numerical  error  including  that 
due  to  artificial  dissipation  which  is  used  to  stabilize  the  solution  procedure. 
Although  there  is  undoubtedly  room  for  further  evaluation  and  development  of  the 
turbulence  model  being  used,  improvement  in  the  numerical  accuracy  with  which  these 
flows  can  be  computed  is  needed  before  this  can  be  undertaken. 
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